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PETSc’s DMPlex interface for unstructured meshes has been extended to support non-conformal meshes. The 
topological construct that DMPlex implements—the CW-complex—is by definition conformal, so representing non- 
conformal meshes in a way that hides complexity requires careful attention to the interface between DMPlex and 
numerical methods such as the finite element method. Our approach—which combines a tree structure for subset- 
superset relationships and a “reference tree” describing the types of non-conformal interfaces—allows finite element 
code written for conformal meshes to extend automatically: in particular, all “hanging-node” constraint calculations 
are handled behind the scenes. We give example code demonstrating the use of this extension, and use it to convert 
forests of quadtrees and forests of octrees from the p4est library to DMPlex meshes. 
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1. INTRODUCTION 


PETSc Balay et al. 2015a; Balay et al. 2015b is an actively developed and widely used library 
for numerical methods in scientific computing, providing parallel data management, structured and 
unstructured meshes, linear and nonlinear algebraic solvers and preconditioners, time integrators 
and optimization algorithms. Many of these methods (such as geometric multigrid and domain 
decomposition linear and nonlinear solvers) can take advantage of the geometric/topological set¬ 
ting of a discretized problem, i.e. mesh information. PETSc’s interface for serving mesh data to 
numerical algorithms is the DM object. PETSc has native DM implementations for several mesh 
formats, and imple mentations tha t wrap external libraries may also be registered, such as DM- 
MOAB for MOAB Tautges et al. 2004 . Because of PETSc’s pointer-to-implementation approach 
to method extensibility, externaFTniplementations may cover only those methods in the DM API 
that are necessary for their target applications. While no DM implementation is privileged above 
others— PETSc-native and external implementations are registered in the same way—the native 
implementations of structured grids (DMDA) and unstructured meshes (DMPlex) have the most 
complete coverage of the DM API, and are developed most actively. Only DMPlex, for example, cur¬ 
rently has complete support for the PetscFE and PetscFV implementations of the finite element 
method and finite volume method. 

Many mesh formats lie between structured grids and unstructured meshes and can broadly be 
described as hierarchical mesh formats. Examples include red-green refinement of triangular meshes, 
quadtree/octree refinement of quadrilateral/hexahedral meshes, and nested Cartesian grids. These 
formats are often implemented in frameworks with data structures that are advantageous for certain 
data access patterns. Patch-based nested grids, for instance, are optimized for fast stencil operations. 
Another example is red-green refined triangular meshes, where the triangles have been ordered by a 
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Sierpinski curve: these meshes can efficiently compute residuals of discontinuous-Galerkin or finite 
volume operators without explicitly determining which cells are adjacent to each-other by pushing 


fluxes onto stacks Bader et al. 2012 


We would like to be able to convert hierarchical mesh formats into the DMPlex unstructured 
format. Our main reason is that the efficiency gains of hierarchical mesh formats typically come at 
the expense of flexibility and generality. DMPlex supports, for instance, arbitrary mesh partitions 
and the extraction of arbitrary subsets of cells (or facets) as submeshes: features which are typically 
missing from hierarchical meshing frameworks. The broad support of DMPlex for PETSc’s DM API 
also makes it an ideal format for testing and comparing numerical methods that call on the DM 
interface. 

What has prevented the conversion of these meshes in the past is that hierarchical meshes are 
often non-conformal meshes: this is true of quadtrees, octrees, nested Cartesian grids, and of some 
hierarchical simplicial meshes as well. CW-complexes—the topologies DMPlex was designed to 
represent—are by definition conformal. Our recent extension of DMPlex has addressed this short¬ 
coming. In this paper we describe how we represent non-conformal hierarchical meshes in DMPlex 
in a way that minimally disturbs the way DMPlex interacts with the other components of PETSc, 
and that requires minimal input from the user. 


2. PRELIMINARIES: CONFORMAL MESHES 

We begin with a brief review of the DMPlex interface and the finite element method in the context 
of conformal meshes. 


2.1. CW-complexes and DMPlex 

The triangulation of a domain D into cells generates a CW-complex (see, e.g., Hatcher 2002 
Chapter 10]). In short, a CW-complex is a partition of a d-dimensional spaces into weil-shapec 
open cells with dimensions between 0 and d, such that the boundary of each n-cell (n > 0) is 
partitioned by finitely many lower-dimensional cells. We call cells of every dimension “points” in 
the complex. 

In a CW-complex, the basic relationship that defines the topology is the map from an n-cell A to 
th e (n — l)-cells on its bou nda ry, which we call the cone of A , cone(A), following the terminology 


Knepley and Karpeev 20091. The closure of the cone map, 


clos(A) := {A} U cone(A) U cone(cone(A)) U ..., 


(1) 


corresponds to closure in D, i.e., clos(A) partitions A. 

The reverse map, taking the n-cell A to its adjacent (n + l)-cells, is called the support map, 
supp(A), and the closure of the support map is called the star of A , star(A). It is important to note 
at this point that for conformal meshes, cones and supports are dual, 

B € clos(A) A e supp(H). (2) 

A CW-complex can be represented by a Hasse diagram for stratified partially-ordered sets: the 
depth of a stratum corresponds to the topological dimension of its points, upward arrows represent 
cone maps, and downward arrows represent support maps. These concepts (cones, supports, strata) 
are at the core of the DMPlex interface, which we illustrate in Figs. [I] to [5] 


2.2. The reference element and element maps 

We assume a Ciarlet reference finite element (K, P(K ), E) (reference cell, space, and dual basis) is 
specified and a domain D is triangulated into a mesh of Nk cells, with cell Ki being the image of K 
under a smooth embedding tpi : K —¥ £1. A finite-dimensional subspace V), of a function space K(f2) 
is then specified as the set of functions v € V such that the pullback ip*v := voipi is in P(K) for each 
ifii. (Discretizations of and iJ dlv (f])-conforming spaces are often pulled back onto reference 


ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY. 







Support for Non-conformal Meshes in PETSc's DMPlex Interface 


A:3 




Fig. 1. A two-triangle mesh and its Hasse diagram. 




Fig. 2. cone(A) / DMPlexGetCone(). 




Fig. 3. clos(A) / DMPlexGetTransitiveClosure(useCone=PETSC_TRUE). 




Fig. 4. supp(£) / DMPlexGetSupport (). 
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Fig. 5. star(<5) / DMPlexGetTransitiveClosure(useCone=PETSC_FALSE). 




Fig. 6. The reference triangle, reference quadrilateral, and the Hasse diagrams of their CW-complexes. 


elements using c ovar iant and contravariant Piola transformations, which have minor implications 
discussed Section 6.1 ) The adjoint of the pullback is the pushforward (<pi*cr)(v) := a(ip*v): it pushes 


E forward onto a set of functionals in V f *, 


Hi := (3) 

To the conventional triplet (K,P(K),'E), we add a reference CW-complex S, which decomposes 
the closure of K (Fig. [Hj) . 

2.3. The finite element method for conformal meshes in DMPlex 

We now list three assumptions that are implicit in most finite element discretizations for confor¬ 
mal meshes, and thus in the data structures and functions that DMPlex uses to implement the 
finite element method. Making these typically implicit assumptions explicit will help to explain the 
extension of DMPlex to non-conformal meshes in Section [3] 

First, we assume that the reference complex S also decomposes the dual basis E, in that the 
shape function associated with each er £ E is supported in the star of a point in S. We formalize 
this as assumption I. 
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I. For each cr*, € E there is a point p £ S such that, if ipk £ P(K) is oVs shape function (crj(ipk) = 
5jk), then supp (ipk) = Ustar(p). 

Assumption I is satisfied by essentially all finite elements: it allows for the definition of compactly 
supported basis functions of Vh- We will refer to the reference functionals associated with p £ S as 
E p , so that E = U pg gE p , and to the pushforward of those functionals under p, as E p := so 

that E i = U pgi jE p . 

Second, we assume that the embeddings of neighboring cells are compatible, in that the traces 
of their approximation spaces “line up” so that H l (VL) functions can be constructed, which we 
formalize as assumption II. 

II. If C ■.= Ki n Kj ^ 0, then ip £ P(p~ l (C)) => tp*ip~*ip £ P(^“ 1 (C)) (where ipj* := (p*) -1 and 
P{X) is the trace space of P(K) on X C K). 

Finally, we assume that the dual bases of adjacent cells are compatible, in that the mappings of 
adjacent cells push functionals forward on top of each other, which we formalize as assumption III. 

III. If p, q £ S and there are adjacent cells Ki and Kj such that ipi(p) = <Pj(q), then there is a 
permutation M such that E p = MEj. 

The permutations typically encode the symmetries of the polytopes in S, e.g., reversal for edges, 
and dihedral symmetries for faces. 

For each vector v £ 14 and each element K t , we need to be able to evaluate E i(v). Given a 
choice of basis IT for V£, each element has a restriction matrix Ri such that Ej(u) = RiW(v). For 
a conformal mesh, assumptions I, II, and III allow for a global nodal basis to be defined for V£ : a 
basis W that is the union of the pushforward dual bases, 

W:= USU^EJ. (4) 

We may also think of IT as being decomposed into the functional associated with points in S , 

W = U seS E s , (5) 

where E s := E p (up to a permutation) if (pi(p) = s. With a global nodal basis, the restriction matrix 
Ri for each cell Ki is a binary matrix, and there is a subset IT, of W such that E* = RiW = IT,. 

Here we see the utility of representing a conformal mesh as a CW-complex S. Given a map 
G : S —> 2 m from each point in S to its set of associated functionals in IT, we can compute IT, as 
G(clos(ifj)), since clos(JG) is the image of the reference complex S under ipi. 

In DMPlex, the map G is represented by a PetscSection. For a typical finite element, the number 
of functionals associated with a point p £ S is a function only of p’s topological dimension, so that G 
can be calculated purely from the sizes of the strata of S': if a DMPlex has been given a finite element 
object (PetscFE), it constructs G automatically, and makes it available by DMGetDefaultGlob- 
alSectionQ. The set clos( A,) can be constructed with DMPlexGetTransitiveClosure(), which 
is used to construct E i(v) on a vector v £ V)* in the function DMPlexVecGetClosure(). This 
function is called within tight, performance critical loops when computing residuals or calculating 
Jacobians (which we illustrate in a prototypical residual evaluation function in Ex. IT]) , so when con¬ 
sidering representations of non-conformal meshes in DMPlex, we chose to avoid those that would 
require modifications at this level of granularity. 

3. NON-CONFORMAL MESHES 

In this section we describe the way non-conformal meshes can now be represented in DMPlex, and 
we describe a general approach to computing with finite elements on these meshes. 
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EXAMPLE 1: finite element residual r = /(v) using DMPlex (dm) 

1 Petsclnt c, cSize, cStart, cEnd; 

2 Vec vLocal, rLocal; 

/* Get vectors for the local representations */ 

3 DMGetLocalVector (dm,&vLocal); 

4 DMGetLocalVector (dm,ftrLocal); 

/* Get the support of v on the local subdomain of this MPI process */ 

5 DMGlobalToLocalBegin(dm, v,INSERT_VALUES,vLocal); 

6 DMGlobalToLocalEnd (dm,v,INSERT .VALUES,vLocal); 

/* Get the range of cell indices (cells are the lowest stratum) */ 

7 DMPlexGetHeightStratum(dm, 0,&cStart,&cEnd); 

8 if (cEnd > cStart) { 

9 PetscScalar *vElem, *rElem; 


/* Get the size of the element dual space £ */ 

10 DMPlexGetVecClosure (dm,NULL,vLocal,cStart,&cSize,NULL); 

/* Get workspace arrays */ 

11 DMPlexGetWorkArray (dm,cSize,PETSC_SCALAR,&vElem); 

12 DMPlexGetWorkArray (dm,cSize,PETSC_SCALAR,&rElem); 

/* Compute the local residual */ 

13 VecSet (rLocal,0.0); 

14 for (c = cStart; c < cEnd; C++) { 

/* Get the restriction of v to cell c */ 

15 DMPlexVecGetClosure (dm,NULL,vLocal,c,&cSize,&vElem); 

/* Compute the element contribution to the residual */ 

/* ... [rElem = f(vElem)] */ 

/* Sum the element residual into the local residual vector */ 

16 DMPlexVecSetClosure (dm,NULL,rLocal,c,rElem,ADD .VALUES) ; 

17 } 

/* Free workspace */ 


18 DMPlexRestoreWorkArray(dm,cSize,PETSC_SCALAR,&vElem); 

19 DMPlexRestoreWorkArray (dm,cSize,PETSC.SCALAR, fcrElem); 

20 } 

/* Sum process contributions into r */ 

21 VecSet (r,0.0); 

22 DMLocalToGlobalBegin (dm,rLocal,ADD.VALUES,r) ; 

23 DMLocalToGlobalEnd (dm,rLocal,ADD.VALUES,r) ; 

/* Free local vectors */ 

24 DMRestoreLocalVector (dm,&vLocal); 

25 DMRestoreLocalVector (dm,ftrLocal) ; 


3.1. Representing non-conformal meshes in DMPlex 

The representation of non-conformal meshes that has been added to DMPlex is limited to hierar¬ 
chical non-conformal meshes. By hierarchical we mean that two mesh points p,q £ S overlap only 
if one is a superset of the other, 


(png/fl)^((pCg)v(gC|))). (6) 

Constructing function spaces on non-conformal meshes that are not hierarchical is more compli¬ 
cated, and not considered here. 

In Fig. [T] we show a simple three-triangle mesh with a non-conformal interface between triangle 
A on one side and triangles B and C on the other, and we have also labeled all of the edges and 
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Fig. 7. A simple non-conformal mesh. 


vertices of the triangles. The question when representing this mesh in DMPlex is how to treat the 
edge c and the edges d, e, and vertex S that overlap it. 

The only way to represent this mesh as a true CW-complex (which cannot have overlapping 
points) is to remove the long edge c (Fig. [ 8 ]), but then the cell A would not be considered a triangle: 
it would be considered a degenerate quadrilateral, with d and e in the cone of A. This is a poor 
format for finite element computations for two reasons. The first is that the support of a basis 
function is no longer correlated with the star operator. A hat function centered at vertex 7 , for 
instance, is non-zero on cell B , but B star( 7 ). The second reason is that the shape of clos(A) is 
not the same as for other triangles, so DMPlexVecGetClosureQ would require special handling to 
restrict a function to A. 

Given these considerations, we include both super-points and sub-points (or “parents” and “chil¬ 
dren” henceforth) in the representations of non-conformal meshes in DMPlex. In addition, we make 
the following extensions to the format, illustrated in Fig. [9] 

(1) We add to the Hasse diagram, which is traversed with cone() and suppQ operations, a sep¬ 
arate tree structure, which is traversed with parentQ and children() operations (in DMPlex, 
DMPlexGetTreeParent() and DMPlexGetTreeChildrenQ). 

(2) We break the duality between cone() and suppQ operations. In particular, the cone of an n-cell 
p includes the “natural” decomposition of its boundary into (n — l)-cells (e.g., the three edges 
of a triangle), while the support of p is the set of (n + l)-cells whose boundaries intersect p. It 
is still the case that q £ conefp) => p £ supp(g), but the converse is not true for non-conformal 
meshes. In the mesh in Fig. [9| for example, A £ supp(d) because dA n d ^ 0, but d qL cone(A), 
because it is not one of the canonical edges of A. The extra support maps are included because: 

It ensures that star(p) covers the support of p’s basis functions, which is important for 
determining the sparsity pattern of finite element matrices (if star(p) D star (q) ^ 0 , then 
there may be non-zeros entries in a finite element matrix for their degrees of freedom). 

It ensures that the support of a facet (a (d — l)-dimensional cell) can be used to identify 
neighboring cells for finite volume and discontinuous Galerkin methods. 
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Fig. 8. 


A true CW-complex representation of Fig. [7] demonstrating that star( 7 ) does not include B. 


7 




Fig. 9. An illustration of the extension of DMPlex for non-conformal meshes. The snaking lines show parent (d). 
parent(e), and parent(5); the bold support arrows do not have matching cone arrows, breaking the duality that is 
present in conformal meshes. 


These extensions pass the minimum bar of not affecting the behavior of DMPlex for conformal 
meshes, but what we really want is for future extensions of DMPlex designed for conformal meshes 
to work automatically for non-conformal meshes as well . Because the P ETSc developers encourage 
contributions from users, including to DMPlex (see |Lange et ah 2015 for a recent example), this 
requires careful attention to the modifications to the DMPlex interface, which we will discuss in 
the next section. 
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Fig. 10. Even though the embeddings and ipj of K r and Kj are curvilinear, <^ . 1 o if,, is component-wise affine, 
so H 1 (Q) -conforming spaces can be constructed from tensor-product polynomials on K. 


3.2. The finite element method for non-conformal meshes 

For a non-conformal mesh, a global nodal basis as defined in the previous section is generally not 
possible: the union of all element functionals, 


:= US u Ef 


pes 


( 7 ) 


will contain linear dependencies. For a hierarchically non-conformal mesh 5, however, it is possible 
to construct a global basis W c that is nearly nodal, by including only the functionals of points that 
have no ancestors, 


W c — i \ Nk i i yV 

1=1 {p € S : parent(<pj(p)) = 0} ' 


( 8 ) 


There is then a constraint matrix 7“ such that W u (v) = I™W c (y) for all v £ Vh- These are 
sometimes referred to as “hanging-node” constraints. In this section, we describe the general method 
for calculating /“. 

We retain assumptions I, II, and III from Section [2] when considering non-conformal 
meshes. Assumption II -that neighboring approximations spaces “line up” for 77 1 (fl)-conforming 
constructions—limits the types of non-conformal interfaces that can occur. Given neighboring cells 


Ki and Kj and p,q £ S such that <fii(p) C then p*<p~* must map P(ipJ 1 (pi(p)) onto P{p). 

For simplicial elements with polynomial spaces, this typically means ipj 1 o ^p i : p —> q is affine; 
for hypercube elements with tensor-product polynomial spaces, pj 1 o ip i must be component-wise 
affine. We note that this is not a requirement that ipi or ipj be affine (Fig. 101. 

When these conditions are met, we can expand each functional in E? in terms of functionals in 


Ej: for each oy £ E p , 


r r)(v) = (ipi*a r )(ip~*ip*v) (9) 

= (v?7*V«oy)(<A^) (10) 

= (vj,^<Pi*^r)(lpsWs(V l jV) ( 11 ) 

= E ( l Pj* l Pi* a r)(' l Ps)^s{v), (12) 

(7 s G Ej 
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Fig. 11. The reference tree T for red-green refinement of triangles. Notice that the reference element S is included 
as a sub-complex of T. 


where ip s is again the shape function of oy. The transferee! functional p^ pi*o r can only be sup¬ 
ported on q, so {p~* Pi*<J r )(ip s ) = 0 if q ^ supp(-!/y). By assumption I, this means that the terms in 
Eq. (12) are non-zero only if ay £ E* for t £ clos(parent(yy (p))), 


(</y*oy)(u) = 


E 


{Pj*Pi*Vr)(lp S )crs(v). 


(13) 


Ts £ Uteclos(parent((^i(p)))^ 


Equation (131 illustrates the two key points needed to compute the constraint matrix If 


(1) If p £ S' and parent(p) ^ 0, then functionals in W u associated with p are linear combinations of 
the functionals associated with points in clos(parent(p)). If any of the points in that set has a 
parent, then we can iteratively apply closo parent to find p's anchor points , whose functionals 
will be in the global basis W c . This lets us compute the sparsity pattern of If. 

(2) The matrix that interpolates to E p from its anchor points’ functionals has entries of the form 
{pJ*Pi*& r ){ips) for ay £ S and shape function ip s £ P(/\). This lets us compute entries in If. 


Without further information, entries in If must be computed by evaluating transfered functionals 
of the form p~fp iir a for any non-conformally adjacent cells K j and Kj. In practice, however, non- 
conformal meshes are usually generated by some predefined set of mesh refinement rules. These 
rules can be encoded in a small non-c onf ormal mesh that we call the reference tree T . For a mesh 
created by red-green refinement (Fig. 11), for example, every non-conformal transfer map pj 1 o pi 
is like mapping one of the edges of the coarse cell (a, b, c) to one of the refined edges (e, /). Due 
to symmetry, all we have to evaluate are the transfered functionals for (pj 1 o pf) ~ (b i—> e) and 
(pj 1 o p^ ~ (6 h-> /), which can then be copied into the correct locations in If. 


4. THE DMPLEXTREE INTERFACE 

Support for non-conformal meshes in DMPlex is available in the latest release of PETSc (v3.6), and 
full documentation can be found online. In this section we introduce the most important components 
of the interface, starting with the highest-level methods that require the least intervention from the 
user, and descending into some of the finer controls available to experts. 

In Section |3.2[ we described how the existence of a predefined refinement pattern, encoded in a 
reference tree T, can enable DMPlex to compute the constraint matrix If more efficiently. The ref¬ 
erence tree is also represent by a DMPlex that is assigned to the target mesh (Ex. [2]). Reference tree 
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Fig. 12. The refinement patterns for the reference trees created by DMPlexCreateDefaultR,eferenceTree(). 
The reference trees themselves also contain the original coarse cells. 


implementations for isotropic refinement on simplices and hypercubes for d = 1, 2, 3 are provided 
by PETSc (Fig.pl. 


EXAMPLE 2: Creating a simplicial mesh dm and setting the reference tree 

1 MPI.Comm comm = PETSC_COMM_WORLD; 

2 PetscBool isSimplicial = PETSCLTRUE; 

3 DM dm, refTree; 

4 PetscSection parentSection; 

5 Petsclnt ^parents, *childIDs; 

/* Create a DMPlex, using, e.g., DMPlexCreateFromDAGO */ 

6 DMPlexCreate... (comm,...,&dm); 

/* Create a reference tree that describes the type of non-conformal interfaces in the 
mesh */ 

7 DMGetDimension(dm,&dim) ; 

8 DMPlexCreateDefaultReferenceTree (comm,dim,isSimplicial,&refTree); 

9 DMPlexSetReferenceTree(dm,refTree) ; 

/* dm retains a reference to refTree, this reference can be destroyed */ 

10 DMDestroy (ftrefTree); 


One can create a conformal DMPlex mesh from just the cone maps (DMPlexCreate¬ 
FromDAGO), and DMPlex can infer the support maps: likewise, one can create a non-conformal 
mesh from just the parent maps. In Ex. [3j we demonstrate setting up the parent maps for the 
simple non-conformal mesh in Fig. [TJ assuming that the red-green refinement in Fig. [TT] is used as 
a reference tree. 

A DM can encompass not only a mesh, but also the fields discretized on it, using a common 
interface for both finite element and finite volume methods. In Ex. [1J we set a standard Lagrange 
Pi (A') finite element on a mesh. 

With the reference tree (DMPlexSetReferenceTree()), the parent maps (DMPlexSet- 
TVeeQ), and the finite element (DMSetField()), PETSc will: 


determine the size of the global vector space (the size of W c in Eq . <©)> 
compute the constraint matrix 7“ from point constraints (Eq. (131), 

apply 7“ when getting the local form of a vector (Ex. I] linepl), so that the local form represent the 
vector evaluated at the unconstrained functionals W u (Eq. (7|), and DMPlexVecGetClosure() 
gets the vector evaluated at the element functionals 10 (Eq. ®o, „ _ 

apply 7“ T when combining local residuals into a global residual (Ex. [l] line 23), 


ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY. 




























A:12 


T. Isaac and M. G. Knepley 


EXAMPLE 3: Setting the parentQ maps and child IDs for dm 

/* The figures use symbols for points, but we have to assign numbers to them. We count 
across each stratum, starting at the bottom. */ 

1 Petsclnt numPoints =16, c=5, d=6, e=7, delta = 14; 

/* c is the parent of each of the children */ 

2 Petsclnt parents[3] = {c, c, c}; 

/* Set numbers for the relevant points in the reference tree as well. */ 

3 Petsclnt bRef = 4, eRef = 7, fRef = 8, deltaRef = 12; 

/* the childIDs are the points in the reference tree to which the children are 

analogous, d is to its parent (c) as eRef is to its parent (bRef), so that is its 
childID. */ 

4 Petsclnt childIDs [3] = {eRef, fRef, deltaRef}; 

5 PetscSection pSec; 

6 MPI_Comm comm = PetscObjectComm((PetscObject)dm); 

7 PetscSectionCreate(comm,&pSec) ; 

8 PetscSectionSetChart (pSec,0.numPoints); 

9 PetscSectionSetDof (pSec,d,1); 

10 PetscSectionSetDof (pSec,e,1); 

11 PetscSectionSetDof (pSec,delta,1); 

12 PetscSectionSetUp(pSec) ; 

13 DMPlexSetTree (dm,pSec,parents,childIDs); 

14 PetscSectionDestroy (&pSec); 


EXAMPLE 4: Adding a finite element to a mesh dm 

/* We are creating a scalar field */ 

1 Petsclnt numComp = 1; 

/* We have a simplicial reference element */ 

2 PetscBool isSimplex = PETSC_TRUE; 

/* The options prefix, for setting options at runtime: e.g., one can change the 

approximation order with ‘-my_fe_petscspace_order 2‘ */ 

3 const char ^prefix = "my_fe_"; 

/* The quadrature order */ 

4 Petsclnt qorder = 1; 

5 Petsclnt dim; 

6 PetscFE fe; 

7 DMGetDimension(dm,&dim) ; 

8 PetscFECreateDefault (dm,dim,numComp,isSimplex.prefix,qorder,&fe); 

9 DMSetField(dm,0, (PetscObject) fe); 

10 PetscFEDestroy (&fe); 


— transform element matrices by the constraints in DMPlexMatSetClosure() to correctly as¬ 
semble a global system matrix from element matrices. 

To create a non-conformal mesh that uses a different refinement pattern than the ones provided 
by PETSc, the user can create a custom reference tree. Any DMPlex that has had the parent maps 
set with DMPlexSetTree () can serve as a reference tree. 

If the user does not provide a finite element, then DMPlex cannot determine for itself the layout 
of the vector space, and the entries in the constraint matrix /“ cannot be calculated automatically. 
If the user specifies the number of degrees of freedom associated with each process-local mesh point 
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(using DMSetDefaultSection()) then the tree data can be used to both compute the size of the 
global vector space and the sparsity pattern of The user can then fill the entries of /“ manually. 

There are also potential uses for intra-mesh constraints between degrees of freedom that do not fit 
into the hierarchical non-conformal framework that is our focus here. The constraint matrix /“ can 
be a PETSc Mat of any specification, and can be added to a DM directly with DMSetDefault- 
Constraints(). These constraints are applied at the conclusion of DMGlobalToLocalEnd(), 
which gets the process-local representation of a vector, and the transpose of these constraints are 
applied at the beginning of DMLocalToGlobalBegin(), when the contributions of all processes 
are summed into a global vector. 

5. VERIFICATION AND EXAMPLE USAGE 

In DMPlex’s example program ex3[j] we include a small verification that DMPlex handles non- 
conformal meshes properly. The example can be run to create simplicial or hypercube meshes with 
non-conformal interfacesQjj To test that the finite element computations are handled correctly, we 
construct a symmetric-gradient Laplacian operator E , 

E(u, v) := / |(Vu + Vu T ) : |(Vw + Vv T ) dx, (14) 

Jn 

and check whether rigid-body motions are in the null-space of E. The rigid-body motions will be 
in the null-space of each element matrix that is computed, but if /“ is incorrect, then they will be 
summed into the global system matrix incorrectly, and it is very unlikely that they will be in be in 
the null-space of the incorrect matrix (Ex. [5]). 


EXAMPLE 5: Testing the correctness of an assembled Jacobian for a non-conformal mesh dm 
(abridged from ex3.c) 

1 MatNullSpace sp; 

2 Vec local; 

3 PetscBool isNullSpace; 

/* This tests that the global system size is determined correctly, and that the 

sparsity pattern for global system matrices is computed correctly */ 

4 DMCreateMatrix(dm,&E) ; 

5 DMGetLocalVector (dm,&local); 

/* This is a finite-element loop within PETSc’s SNES library that assembles the 
Jacobian matrices of nonlinear equations: the vector local is needed as a dummy 
argument to represent the current "solution" used to evaluate the Jacobian, which in 
this case is the linear operator in Eq. ( | 14| > . The variational form needed to compute 
each element’s matrix has already been attached to dm. */ 

6 DMPlexSNESComputeJacobianFEM (dm,local,E,E,NULL); 

7 DMPlexCreateRigidBody(dm,&sp) ; 

8 MatNullSpaceTest(sp,E,&isNullSpace) ; 


1 src/dm/impls/plex/exEmiples/tests/ex3. c: run make ex3 in that directory to build the example. 

2 The initial intent of our work is merely to allow non-conformal meshes to be represented in DMPlex, not to 
implement a stand-alone adaptive mesh refinement interface. To test the DMPlexTree interface without relying on 
external libraries, however, we have written DMPlexTreeRefineCell(), which hierarchically refines a single cell of 
a conformal mesh. 

3 To visualize the non-conformal meshes used, go to the src/dm/impls/plex/examples/tests/ directory of the 
PETSc source and run make ex3; ./ex3 -tree -simplex B -dim D -dm.view vtk:nonconf_BJ).vtk:ASCII.VTK for 
B S (0,1} and D S {2,3}. 
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Fig. 13. A Mobius strip mesh, generated in p4est and converted to DMPlex in the example program p4est_test_plex. 

Example usage outside of PETSc can be found in the p4est library for parallel adaptive mesh 


refinement Burstedde et al. 2011), which implements the forest-of-quadtrees and forest-of-octrees 


paradigms in 2D and 3D. This library is meant to provide data structures only, and comes with no 
built-in solver or finite element f ramework . By c onverting the p4est format into DMPlex (building 
on the methods described in |Isaac et al. 2015) for efficiently converting p4est’s native format to 
adjacency-based formats like DMPlex), we make PETSc’s numerical methods more readily avail¬ 
able to p4est users. Example pr ograms that perform this conversion are distributed with p4est as 
p4est_test_plex (2D) (Fig. 13) and p8est_test_plex (3D). The repository of the p4est library]^] 
has a “petsc” branch that is compatible with PETSc 3.6. [^] 

6. OTHER DISCRETIZATIONS 

The focus of this work has been if 1 (0)-conforming finite elements. We briefly discuss the way our 
approach to non-conformal meshes in DMPlex affects other discretizations. 

6.1. H curl (Q)- and 4f d '“(f2)-conforming finite elements 

PetscFE does not currently implement the covariant and contravariant Piola transforms that are 
commonly used by H curl (Q)- and if dlv (D)-conforming finite elements, but these methods can still 


be formulated via pullback onto reference elements Rognes et al. 2009 , so future PetscFE im¬ 
plementations of these finite elements are a possibility. The discussion of conformal meshes and 
non-conformal meshes in this work is still valid for these finite elements, with two small modifica¬ 
tions: 


The pullback operations are defined to be 
<P*v '■= Vtpjvotpi 
ip*v := |detV(/?i|V</?iU o 


[covariant, 7T curl (fl)], 
[contravariant, iT dlv (S 2 )]. 


(15) 

(16) 


The trace space P{p) for a point p £ S involves not only restricting the function space P{K) 
to the point, but also restricting to the tangential component (iT curl ) or the normal component 
( H div ). 


4 https: //bitbucket.org/cburstedde/p4est / 

5 To build these examples, run ./configure —with-petsc=$PETSC_DIR and make test/p4est_test_plex 
test/p8est_test_plex. To view the DMPlex meshes created in these tests, run the examples with the flag 
-dm_view vtk:p4est_petsc.vtk:ASCIIJVTK. 
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The core operation to compute continuity constraints is the transfer of a functional from one 
element A'.j to its neighbor Kj and evaluation on a shape function, (¥y t 1 </?i* (T )(V’) = cr(ip*(pj*ip). 

After one verifies that ViPj* = (rj 1 ° Vi)* f° r both pullbacks above, then it must be true that 
if each child-to-parent map ipj 1 o ^p i is represented in the reference tree T, then it can be used to 
compute the entries in the constraint matrix 

6.2. The finite volume method 

Finite volume methods do not promote the encapsulation of complexity as well as finite element 
methods. We have formulated our non-conformal mesh extension for finite elements such that, in 
a typical finite element loop (Ex. fp, the operations performed on each cell in the loop do not 
depend on whether or not any of the points in the cell’s closure is a child or a parent. In a cell¬ 
centric approach to the finite volume method, the act of reconstructing centroid values requires 
determining the neighbors of a cell, which becomes more complex when multiple cells may be on 
the opposite side of a face (in Fig. e.g., both cells B and C are opposite cell A across edge 
c). While we are currently incorporating non-conformal meshes into the finite volume method as 
implemented by PetscFV (we expect to finish while this manuscript is in review), the result is likely 
to be more fragile to user extensions. 

One particular aspect that will be counterintuitive to users who worked with finite volume meth¬ 
ods on conformal meshes is that a facet can have more than two cells in its support. One often 
finds in finite volume code constructs of the form “neighbor = (supp[0] == me) ? supp[l] : 
supp [0],” which are no longer valid. One also has to avoid double-counting fluxes, i.e., computing 
fluxes on both a parent facet and its children. 

Because many unstructured finite volume methods do not care about the shape of cells (i.e., 
whether they are triangles or quadrilaterals), these issues can be avoided by encoding non-conformal 
meshes as conformal (though degenerate) ones, as in Fig. [8] In a multiphysics setting, where a finite 
volume field and a finite element field are involved in a larger system of equations, this approach is 
not possible. 

7. DISCUSSION 

We have presented an extension to PETSc’s DMPlex interface for unstructured meshes so that 
it can now represent hierarchical non-conformal meshes. Our extension leaves the interface for 
conformal meshes the same, but adds a tree structure to encode the hierarchy of subsets (children) 
and supersets (parents). We have shown how, for a wide class of finite elements, by combining this 
hierarchical information with a reference tree that describes the types of non-conformal interfaces 
that appear in a mesh, the extra complexity of non-conformal meshes can be hidden, allowing finite 
element code written for conformal meshes to be applied to them. This extension can already be 
used to convert p4est forest-of-quadtrees and forest-of-octrees meshes to DMPlex. Work is underway 
to bring support for the finite volume method up to the level of the finite element method, and 
future work on the discontinuous Galerkin method is planned. 
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